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Abstract 

The contour integration technique applied to calculate the optical conductivity 
tensor at finite temperatures in the case of layered systems within the framework 
of the spin-polarized relativistic screened Korringa-Kohn-Rostoker band structure 
method is improved from the computational point of view by applying the Gauss- 
Konrod quadrature for the integrals along the different parts of the contour and 
by designing a cumulative special points scheme for two-dimensional Brillouin zone 
integrals corresponding to cubic systems. 



1 Introduction 



Nowadays magneto-optical effects are widely used to probe the magnetic properties of 
various systems |2], [|. F° r a theoretical description of these effects, one needs to 
calculate the optical conductivity tensor in a parameter-free manner. Recently, two of 
the authors have proposed a new, contour integration technique to calculate the optical 
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conductivity tensor for surface layered systems J|. The theoretical framework of the 
present paper is based on this technique. Therefore, in the following the basic concepts 
of this method are only briefly reviewed. 

The starting point of the contour integration technique is the expression for the optical 
conductivity tensor 

as given by Luttinger ||, where u denotes the photon frequency and 5 a finite life- 
time broadening, respectively. The latter accounts for those scattering processes, which 
are not incorporated in a standard band structure calculation, but are present at finite 
temperatures. 

The temperature T enters the expression for the zero wavenumber current-current corre- 
lation function parametrically || 
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via the Fermi-Dirac distribution function /(e), with e m and e n being the eigenvalues of the 
one-electron Hamiltonian corresponding to states labeled by m, n and the J% m (// = x, y) 
the current matrices. In Ref. Q it was shown, that Eq. (Q) can be evaluated by means 
of a contour integration using complex energy values z. Within this technique, a^uo) is 
decomposed into a contour path contribution ajS)(uj) and a contribution, affl(u), arising 
from the Matsubara poles a^P(u>), such that 



^(fj) = aff( U ) + a^( U ). (3) 

As shown in Fig. [I], (u), in turn, consists of the contributions from the contour in the 
upper and lower semi-plane as given by Eqs. (24) and (25) in Ref. One contribution 
to a$y{oj) comes from the n 2 Matsubara poles near and on both sides of the real axis, 
and an other one from the n\ poles situated exclusively in the upper semi-plane, see Eq. 
(26) in Ref. Q. (Note that according to Ref. |f, m = Nt - N 2 and n 2 = N 2 .) Each of 
these contributions (altogether four) is expressed in terms of 

**) = Tr [ JfiG ( z ^ rG (to)\ , (4) 

where J M and G(z{) denote current operators and resolvents, respectively. Application of 
this contour integration technique to compute £^(0;) for ordered or disordered (within 
the framework of the single site Coherent Potential Approximation) layered systems is 
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therefore straightforward [0. Originally, the quantities d^iz^z^) were introduced to fa- 
cilitate the computation of the dc conductivity based on the Kubo-Greenwood formalism 
|| H in case of substitutional^ disordered bulk systems JT0| . 



In the present paper, the Green functions G(z) entering Eq. (^) are calculated by means 
of the spin-polarized relativistic screened Korringa-Kohn-Rostoker (SKKR) method for 
layered systems 0, [12|, [13|] and the matrix elements J^ n (/i = x, y) using the relativistic 
current operator 0. Because of the non-vanishing imaginary part of the complex energies 
z (see Fig. [[]), also the so-called irregular solutions of the Dirac equation have to be 
considered 

To start a computation of the optical conductivity tensor E /UV (o;) for a given frequency uj 
and temperature T, the self-consistently determined potentials of the investigated layered 
system must be known. This means that the bottom energy and the Fermi level ep are 
also known. Thus min(Re z) is fixed at and max(Re z) at e-p + mk^T (m G N, k-Q the 
Boltzmann constant). Due to the fast decay of the Fermi-Dirac distribution function, at a 
given T the computed optical conductivity tensor depends only slightly on the used value 
of m. The results given in the present paper were obtained with m — 8. As can be seen 
from Fig. [I], the value of Im z is 8\ in the upper and 5 2 in the lower semi-plane along the 
contour part parallel to the real axis. Taking, 8j = 2nj5T (j=l,2), with 5t = irk^T, it is 
ensured that the paths parallel to the real axis fall in-between two successive Matsubara 
poles g. 

In order to make use of the symmetry properties of the d , lul {z\^ z 2 ) in Eq. (|j), for the life- 
time broadening a value of 5 = 2# 2 has been chosen. An other advantage of this set-up 
is that for example at T = 300 K and 5 ~ 0.05 Ryd, one needs only n 2 = 2 Matsubara 
poles. 

The computed value of the optical conductivity Ti^iuS) depends strongly on the number of 
complex energy points used on different parts of the contour in both semi-planes: quarter 
circle, parallel to the real axis and in the vicinity of the Fermi level. Furthermore, the 
number of fc-points used to calculate the scattering path operator and o^iz ±hw + i8, z) 
at a given energy z, (see Eqs. (53) and (54) from Ref. 0), is of crucial importance. 

The aim of the present paper, is to describe efficient numerical methods to control the 
accuracy of the contour and fc-space integrations without increasing the computational 
effort: in Section |2] Konrod's extension to the Gauss quadrature as a proper numerical 
method is discussed; in Section || a new, cumulative special points method is presented 
to compute two-dimensional /c-space integrals with arbitrary high precision. Finally, 
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the independence of on the contour path in the upper semi-plane, i.e. on the n\ 

Matsubara poles, is shown separately in Section £|. 



2 Contour integration by means of the Konrod— Legendre 
rule 

The n-point Gauss-Legendre integration rule [|14j] can be used directly to compute a^J (uj) 
in Eq. by transforming the nodes Xk and their weights Wk (k = 1, . . . , n) according 
to the different contour parts (Fig. |l|). The nodes are roots of the Legendre polynomials 
P n (x) corresponding to the weighting function w(x) = 1, x G (—1,1) and are usually 



computed using the Newton method flip] . For the weights Wk, one needs also to compute 
the derivatives of the P n (x) with respect to the argument P^Xk). 

One way to estimate the accuracy is to repeat the above procedure for n + 1 points 
and compare the results. This requires the evaluation of the integrand for all the newly 
generated n + 1 complex energy points, because P n+ i(x) has no common roots with P n (x) 
14], which of course is computationally very inefficient. 

In 1965 Konrod has proposed a method to overcome the above difficulty by demonstrating 
that one can create a set of 2n + 1 nodes including all the nodes of an n-point Gauss 



quadrature [fTE ], Furthermore, he also showed that each of the additional n + 1 nodes falls 



in-between two nodes of the n-point Gauss quadrature. Thus, once the integrand f(x) 
is evaluated for each of the 2n + 1 nodes, Xk, both the Gauss (denoted by Q n ) 

n 

Snf = ^ f(^k)W 2 k (5) 
k=l 

and the Konrod sum (denoted by K,2n+i) 

2n+l 

JC 2 n+lf = f(Xk)w k , (6) 

k=l 

are available. Usually, the nodes Xk and weights, Wk and w 2k , are obtained from the 
spectral factorization of the associated Jacob i-Konrod (Jacobi-Gauss) matrix [I7|. The 



Jacobi-Gauss matrix is formed easily knowing the recursion coefficients of those monic, 
orthogonal polynomials (say Legendre polynomials), which correspond to the weighting 
function (in our case the identity) [14]]. But these coefficients fill only partially the Jacobi- 
Konrod matrix. Hence, to form the Jacobi-Konrod matrix is not a trivial task. 



4 



Laurie |18] was the first, who in 1997 developed an algorithm, based on the mixed moments 
method in order to generate the Jacobi-Konrod matrix for even n. Recently, Calvetti et 



al. [[17]] extended Laurie's ideas to odd n using a divide-and-conquer scheme. 

We have implemented Laurie's algorithm [ll| to calculate the 2n + l recursion coefficients 
entering the Jacobi-Konrod matrix. The initialization requires 3n/2 (n even) elements of 
the Jacobi-Gauss matrix derived from a common set of monic, orthogonal polynomials. 
These are obtained using the algorithm described in Ref. [pTQf] . The same Ref. |l!|] is then 
used to generate the nodes and weights for the Konrod rule and the weights for the Gauss 
rule with machine dependent accuracy, say 10~ 15 . In contrast to the implementation of 
the Konrod rule described in Ref. pOfl , the present one works for arbitrary even n, i.e. Q n f 
can be compared with A^n+i/ not only for some particular values of n. The final result of 
the quadrature is then that obtained by means of the Konrod rule. Since ajS,\uj) is not 
a scalar quantity for each pair of (jj, is), QnO^iuj) = affliu) is compared separately with 
K>2n+i<J iiv{^) = <j^u +1 \uj)- The integrand along a particular part of the contour, see Fig. 
[I], is said to be converged if the maximum difference 

max | a^ +1) M - (w) | < e z (//, v = x, y) , (7) 
where e z is an arbitrary small number, e.g. machine accuracy (10~ 15 ). 

In Fig. | a^u +1 \uo) — affl(uj) | is plotted for fko = 0.05 Ryd and T = 300 K versus the 
number of complex energy points n\_ used parallel to the real axis in the lower semi- 
plane. For all other parts of the contour, not shown here, the corresponding quantities 
| <J^v +1 \oj) — affi(uj) |, show an almost linear dependence on n. The rather complicated 
shape of the data displayed in Fig. [2] has a physical reason: in the lower semi-plane, near 
(62 = 0.024 Ryd) and parallel to the real axis mostly the joint density of states with all 
its singularities is mapped. 

In the present paper, the layered system used for test calculations is a mono-layer of 
Co on the top of fcc-Pt(100), i.e. below the Co surface layer there are three Pt buffer 
layers followed by Pt bulk [^TJ. The band bottom energy (et,) and the Fermi level (ep) of 
the substrate corresponds to —1 and —0.039 Ryd, respectively. The optical conductivity 
calculation was carried out using rii = 18 and n 2 = 2 Matsubara poles. The number of 
/c-points within the surface Brillouin zone was 16 (further discussions on this point are 
made in Section |^). Analyzing the values obtained, it can be concluded that for the case 
chosen, the following set of numerical parameters 

n™ c = 4 77J4. = 8 n £ z F + = 6 (upper semi-plane) 
n circ _ 4 tv\_ = 88 n £ z _ = 8 (lower semi-plane) 
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yields a maximum differences | a^J 1 (u) — aj^J{uj) | < 10 6 a.u. for each part of the 
contour (dotted line in Fig. |2|). Notice that e z = 10 -6 a.u. is hundred times smaller than 
the smallest contour part contribution to (cu). 



3 Cumulative special points method for two— dimensional 
lattices 

In the present paper the special points method (SPM) has been used. The reason for this 
is twofold. It can be shown that the SPM in fact is a Gauss quadrature |[22|| . Therefore, 
its application for a computation of E^cu) guarantees that all integrations involved are 
performed by means of the same quadrature method. Furthermore, as demonstrated 
below, it is possible to use the SPM cumulatively, which in turn facilitates to monitor the 
accuracy of the fc-space integration. 

The integral S Qj , of a function f(k) over and normalized to the surface Brillouin zone 
(SBZ) is approximated in the SPM similar to Eq. (|5|) by 

= ( 8 ) 

i=i 

where n^ denotes the number of special points in the irreducible part of the surface 
Brillouin zone (ISBZ), and the weights Wj have to fulfill the requirement: 

n k 

5>i = l. (9) 

i=i 

These special points kj are defined by the following condition 

^TivU^-Ho, (io) 

namely, in terms of a homogeneous system of linear equations in symmetrized plane 
waves A m (k) [[23|, |25[, which form a set of real, orthogonal, translationally and (point- 
symmetry group) rotationally invariant functions . 

Although there are several methods known in the literature to solve Eq. (|H]) for three- 
P5] , |27| ] and two-dimensional |2S| Brillouin zones, in the following we adopt the scheme 
proposed by Hama and Watanabe They have shown that the set of fc-points 

2 

kj = ^2k ja b a , (11) 
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with 



Ja 



ja ~ 1 



2 ' Jc 



1. 



a 



1,2 



(12) 



are solutions of Eq. (flTf) , i.e. special points, which minimize the remainder in Eq. (||). 
Hence the special points form an uniform, periodic mesh with respect to the edges b Q 
(a = 1, 2) of the reciprocal unit cell |2"7|] , but they are not uniquely defined because of the 
arbitrariness [|3D] of the parameter a a (a = 1, 2) in Eq. (TT2]). 

The extension of the SPM proposed in the present paper exploits the arbitrariness of 
a a and is based on the observation that successively denser fc-meshes, including all the 
Appoints of the previous meshes, can be created, if the parameter a Q in Eq. (|T2| ) does not 
depend on n a . Consider a two times denser mesh 



2n n a 2 



1. 



2n a \ a 



1,2 



(13) 



than that in Eq. (|12|)- This new mesh includes |3T] all the (former) kj— points and has 
additional points in-between, because 



for j c 



, n a and 



--k 



3a ■ 



Ja 1 



Akj 



Akj 



if L 



if L 



^3a 
2 Ja 



(14) 



Ak ia = k 



ia + l 



k: 



It should be noted that the validity of the above statements does not depend on the 
dimensionality of the Brillouin zone. 

In our, cumulative SPM, we use origin centered fc-meshes, i.e. Eq. (O) for a Q = 1/2 and 
the same number of divisions n a = n in each direction a. Since our interest in evaluating 
S M ^(co>) is mainly restricted to cubic layered systems, in Table [I] all details regarding 
the origin centered fc-meshes for primitive rectangular, square and hexagonal lattices are 
listed. 

As long as the magnetization is perpendicular to the surface, the irreducible part of the 
SBZ is identical to the paramagnetic one. (This situation pertains to the present paper.) 
The construction of the paramagnetic ISBZ follows closely the one, introduced years ago 



by Cunningham [28]. However, in the case of the hexagonal lattice, a two times bigger 
ISBZ was taken, obtained by rotating clockwise by 60° his ISBZ and subsequent 
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lattice 


kj e fiisBZ, if 






primitive 
rectangular 


n 

1 < J* < - + 1 

(a = 1,2) 


4-2 (5^1 + (5 i2 i) + ^i^i 


(n + 1) 2 


square 


i < h < h < \ + 1 


4 + 4(1 -5 jlja ) 

(1 - dfcl) - 3^1^! 


(n + 1) 2 


hexagonal 


1 < J'l < H < \ + 1 
Ji + 2j 2 < ri + 3 


6-(25 Jll + 3)5 i2l 


1 - 9[to] ([to] + 1) 
+3 Q + l) (f + 2 M ) 



Table 1: Origin centered two-dimensional fc-meshes for even n in the case of a primitive 
rectangular, square and hexagonal lattice. n^{n) is the number of points in the SBZ 
obtained by dividing the edges of the two-dimensional reciprocal unit cell into n equal 
parts, ([m] denotes the integer part of (n + 2)/6 and Sij is the Kronecker symbol.) 

mirroring along the k x axis. The so obtained ISBZ and /c-meshes are in accordance with 



those of Hama and Watanabe used for the three-dimensional hexagonal lattice. 

The weights Wj in Table [l] were deduced using the elements of the corresponding point- 
symmetry groups, i.e. , C 2v (primitive rectangular), C 4v (square), and C 3v (hexagonal), 
respectively [J32[| . They are normalized to the total number of equivalent /c-points in the 
SBZ (last column of Table H) and fulfill the condition in (^). It should be noted that all 
formulae in Table [I] are valid only for even n. 

When the cumulative SPM is used, 



(15) 



new fc-points are added to a previous fc-mesh (i > 1) and their contribution to the SBZ 
integral to be evaluated is labelled by AS ni f. If no previous mesh exists, n%(iii-i) points 
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are created according to Eqs. ( |TT1) and (|12|) leading to the following normalized sum 



n n (m-i) 

^ g (16) 

As a starting mesh n = 2m (m arbitrary even number) can be used. The subsequently 
created meshes then correspond to rii = 2*n (i = 1,2,3,...) divisions along each vector 
b a . Eq. (|T6|) with An^rii) newly created fc-points can also be used to compute AS ni f. 
Proceeding in this manner, one obtains a recursion relation of type 

SnJ = — 7-r [Sm-J ■ nk( n i-i) + A5 «J ' An k( n i)] • ( 17 ) 

n k\. n i) 

The fc-points to be added to a previous mesh are selected in terms of Eq. (p^D, by imposing 
that in Table [I] ji and j'2 cannot be simultaneously odd. It should be noted that expressions 
to evaluate An^(nj) directly can be also deduced from the n^n) listed in Table [3]. Eq. 
fll?]) is then repeated until the absolute difference between S n J and S ni _ x f is smaller than 



a desired accuracy or an allowed maximum number of /c-points nl max ' ) is reached. In 
particular, for ^^(u), this means that 

max I (T<p\u) - \ < e % (ji, u = x,y). (18) 

is imposed for each complex energy on the contour and Matsubara pole. 

It should be noted that Eq. applies to the full SBZ, whereas Eq. (|]) refers to an 
irreducible wedge of the SBZ 0, 

An application of the cumulative SPM is shown in Fig. || For these calculations (Huu 
= 0.05 Ryd and T=300 K), the same layered system is considered as in Section |[ The 
results obtained with a starting fc-mesh consisting on 15 /c-points in ISBZ (rij_i = no = 8) 
is taken as reference. In Fig. |3] these data are compared with those obtained using 45 
A;-points in ISBZ (n, = 2n = 16). For the contour integrations an accuracy e z = 10~ 4 
a.u. was achieved on each contour part (see Fig. [I] and Section |2|). This means that even 
the minimum of ^^(zi, z 2 ) has a last digit exactly computed. 

As can be seen in Fig. |3|, a common precision of = 10~ 4 a.u. can be achieved easily 
with a one-step cumulative SPM for all parts of the contour and for the Matsubara poles 
situated in the upper semi-plane far off from the real axis. Obviously, for the Matsubara 
poles near the real axis more /c-points are needed in order to achieve the same accuracy 
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4 Contour path independence 



For the layered system described in Section | in Fig. [| the optical conductivity E M1/ (u;) 
[a. u.] for fijjj = 0.05 Ryd and T = 300 K is shown as function of the Matsubara poles 
rii used in the upper semi-plane far off from the real axis. (712 = 2 Matsubara poles were 
used near the real axis in both semi-planes.) The convergence criteria (|7|) and fll8D were 
satisfied for e z = e^= 10~ 4 [a.u.]. 

As can be seen from Fig. [|, the contribution coming from the contour (a;) and from the 
Matsubara poles Y$,(uj), see also Eqs. (24) and (26) in Ref. [f|], in the upper semi-plane, 
respectively, depends remarkably on n\. However, their sum does not really depend on 
ni, i.e. 

E^M + E^H = const. ( M ,^ = x,y) 

with an accuracy of n\E^ ~ 10~ 3 [a.u.]. Hence an evaluation of Y^^iuS) does not depend 
on the form of the contour in the upper semi-plane. 

5 Summary 

The computational scheme for the optical conductivity tensor for layered systems has been 
improved numerically. For the contour integration, the Konrod-Legendre rule was applied, 
showing that any desired accuracy e z can be achieved (in comparison with the Gauss- 
Legendre rule). In the case of the A;-space integration, a cumulative special points scheme 
was developed for two-dimensional lattices. This method permits one to perform for each 
complex energy z the fc-space integration iteratively, evaluating the integrand only for 
those fc-points added to a previous mesh. The thus controlled z— and /c-convergence, 
provides independence from the form of the contour in the upper semi-plane with a 
predictable accuracy. 

It should be noted that the described numerical procedures can be used also for other 
approaches or to calculate other physical properties. For example, the cumulative spe- 
cial points method provides an excellent tool to check the fc-convergence of the band 
energy part of the magnetic anisotropy energy. The numerical efficiency in calculating 
the transport properties can be improved in a similar way. 
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Figure 1: Contour path and Matsubara poles used to calculate the zero wavenumber 
current-current correlation function. (n z ± denotes the complex energy points used on 
different parts of the contour. ri\ and n 2 are Matsubara poles, eb is the band bottom 
energy and e-p the Fermi level, respectively.) 
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Figure 2: Absolute difference between a^y +1 \uo) and a$(u) (atomic units) versus the 
total number of complex energy points used parallel to the real axis in the lower semi- 
plane for fku = 0.05 Ryd and T = 300 K. The real part in the difference of the regular 
(left panels) and irregular (right panels) contribution is given by open symbols (o, □) and 
the imaginary part by corresponding filled symbols (•, ■), respectively. The dotted line 
marks a value of 10~ 6 in atomic units.) 
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Figure 3: Maximum absolute difference between a^J (u) and in atomic units 

for a complex energy points N z on a part of the contour or for a Matsubara pole N±, 
N 2 (n = 8, huj = 0.05 Ryd and T = 300 K). Open (filled) symbols are used for the 
contour in the upper (lower) semi-plane. Circles refer to the quarter circles, squares to 
the parts parallel to the real axis and diamonds in the vicinity of the Fermi. A denotes 
the Matsubara poles far off from the real axis in the upper semi-plane and V the others. 
The dotted line marks the value of 10~ 4 in atomic units. 
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Figure 4: Real (open symbols) and imaginary part (filled symbols) of H^u), respec- 
tively, in atomic units for ftw = 0.05 Ryd and T = 300 K against the Matsubara poles 
n\ in the upper semi-plane. (Triangles up are used for the contribution coming from the 
contour and triangles down for that coming from the Matsubara poles. Their sum is given 
by diamond. The convergence criteria (0) and ([18]) for e z = = 10~ 4 a.u. were satisfied 
for each n\ value considered.) 
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